Interlayer Registry Index of Layered Transition Metal Dichalcogenides

Inspired by the fascinating electronic properties of twisted transition metal dichalcogenides, we extend the registry index approach to quantify the interlayer commensurability of homogeneous and heterogeneous interfaces of MoS2, WS2, MoSe2, and WSe2. The developed geometric measure provides quantitative information about their sliding energy landscape with vast mechanical and tribological implications. Furthermore, the registry index is highly suitable for characterizing surface reconstruction in twisted transition metal dichalcogenide interfaces that dictates their intricate electronic and ferroelectric properties. The simple and intuitive nature of the registry index marks it as a powerful computational tool for studying the fascinating physical phenomena demonstrated by these materials.


Density functional theory computational details
The total energy of the considered transition metal dichalcogenide (TMD) bilayers at various stacking modes was evaluated via density functional theory calculations using the plane-wave pseudopotential Quantum Espresso package. 1 The Perdew-Burke-Ernzerhof (PBE) generalized gradient exchange correlation functional approximation augmented by the Grimme-D3 dispersion correction with the Becke and Johnson (BJ) damping 2 and the PSlibrary 3 projector augmented wave (PAW) pseudopotentials were employed. This approach was shown to capture well the interlayer interactions in TMDs. 4 A plane wave energy cutoff of 60 Ry and a k-mesh of 12 × 12 × 1 using the Monkhorst-pack scheme 5 were used for the calculation of the bilayer model. A sufficiently large vacuum size of 10 nm along the normal direction was set to avoid interactions between adjacent bilayer images. For systems exhibiting a finite electrostatic potential difference between the two layers, the dipole correction was applied.
Starting from the relaxed AA' stacking mode of the antiparallel configurations, two types of calculations were performed: (i) rigid shift calculations, where the layers are laterally shifted rigidly with respect to each other along the armchair in steps of 0.1 Å and the total energy is calculated at each step; (ii) vertically flexible shift calculations, where at each step the vertical (z) coordinates of all atoms are allowed to relax prior to the total energy evaluation. The convergence threshold on forces for structural optimization was set to 1 × 10 −4 Ry bohr and the relaxation was performed with a step size of 0.2 Å. These two simulation procedures were repeated for the parallel configuration, where the initial AB stacked bilayer was constructed from the relaxed AA' configuration via a 60° rotation of the top layer. Table S1 summarizes the structural parameters of the TMD homobilayer of MoS2, MoSe2, WS2, and WSe2, and the corresponding six heterobilayer models. Both the lattice constants ( ) and the interlayer distances (ℎ 0 ) of the structurally optimized AA' stacking mode are provided. The interlayer distance is defined as the vertical distance between neighboring chalcogen atom in the adjacent layers. For the heterojunctions, a strained supercell was constructed with a lattice constant taken as the average of the lattice constants of the two isolated monolayers, relaxed at the same level of theory. For several high-symmetry stacking modes of the WSe2 bilayer we performed reference calculation replacing the PBE functional approximation with the Heyd-Scuseria-Ernzerhof (HSE) 6-9 screened-exchange hybrid density functional. This functional, which is more computationally demanding, is known to be successful in describing the electronic properties of homogeneous and heterogeneous layered material structures. [10][11][12][13][14][15][16] As shown in Figure S1, the difference between the PBE and HSE calculated sliding energy barriers is within 0.6 meV/atom, marking the suitability of the PBE reference calculations for the and parameterization. Naturally, given a full reference HSE based dataset, one can readily reparametrize the / . This, however, is not expected to modify the qualitative nature of our results and conclusions.   interlayer configuration (using the same grid as in the rigid shift case) was scanned in intervals of 0.01 Å −1 , and the was evaluated using the optimal radii set. The value of yielding the minimal average deviation from the vertically flexible DFT reference data was chosen as the optimal value. Table S4 lists the average and standard deviation differences between the reference DFT values and the results obtained using the optimal parameter sets.   In Figures          In the calculations presented above and in the main text we used separate parameterizations for the AP and P interlayer orientations. Figure S10 demonstrates that one can obtain quite good agreement between the DFT reference data and the results (with some minor deviations near = 1.2 ) for rigid shifts, even with a single parameter set. Somewhat larger deviations, however, appear for vertically flexible shifts, possibly due to variations in the interfacial electron density between the P and AP configurations. These, in turn, require different effective atomic radii to represent the corresponding Pauli repulsions.

Unnormalized sliding energy curves
As shown in SI section 2 above, the sliding energy profiles of all bilayers considered show very similar features and nearly overlap when presented in normalized form. For completeness, in Figures S11 and S12 we present the reference DFT sliding energy curves without normalization for the various homogeneous and heterogeneous bilayer systems considered in this work, respectively.

Structural relaxation
In the main text, we presented a local registry index ( ) analysis of surface reconstruction in twisted MoS2 bilayers. To this end, supercells consisting of 7 full moiré patterns of twisted MoS2 bilayers with = 0.25° for the antiparallel orientation and = 0.65° for the parallel orientation, were generated using the LAMMPS package. 17 Periodic boundary condition were applied in the lateral directions. In the vertical direction a sufficiently large vacuum size of 10 nm was applied to avoid spurious interactions between adjacent bilayer images. The intralayer interaction was described via the second-generation reactive empirical bond order (REBO) potential, 18 and the interlayer interaction was calculated by a dedicated registry dependent interlayer potential (ILP). 19 The Fire algorithm with a force tolerance of 10 −6 /Å was used to relax the structure (including the box size) keeping the lower sulfur atomic sublayer of the bottom MoS2 layer fixed to mimic a rigid substrate.
The maps presented in Figure 4 of the main text conform well with the lattice domains found in experiments. 20 Specifically, when plotting a one-dimensional profile cross section along a domain wall (see Figure S13) a typical width of ~ 4 nm (~5 nm) is found for the parallel (antiparallel) orientation in good agreement with measured values. 20 Note that the corrugation calculated for the ILP 19 relaxed MoS2 bilayer is somewhat smaller than the sliding profile presented in Fig. S2b above (e.g. ~0.27 versus ~0.4 for the AA stacking mode). This is mainly attributed to the fact that the ILP was parametrized against non-local many-body dispersion corrected Heyd-Scuseria-Ernzerhof calculations, 19 which produce a somewhat larger interlayer distance.
For example, for the AB stacked MoS2 bilayer, the HSE+MBD-NL interlayer distance is ~0.2 Å larger than that obtained using PBE-D3. S13 Figure S13. One-dimensional profile cross section across domain walls in MoS2 bilayers twisted by = 0.25° for the antiparallel (blue) and = 0.65° for the parallel (red) orientations. S14

Comparison to previous MoS2 registry index definition
A GRI for MoS2 was previously defined and parametrized against a local density approximation DFT sliding energy landscape obtained under external pressures of 500 MPa and 15 GPa. [21][22] This parameterization, which considered only the antiparallel interlayer configuration, provided good agreement with the reference data without accounting for Mo-Mo interactions. For completeness, we compare the GRI profile obtained using the previous definition and parameterization with our new DFT reference data for vertically flexible shifts at the antiparallel configuration. The results, appearing in Figure S14, show good agreement between the previous GRI definition and the new reference DFT data with the former presenting more sharp features due to the use of sharp circles (rather than smooth Gaussian) to represent effective atomic volumes. We note however, that the previous definition is insufficient to provide a good description of the sliding energy profile of the parallel configuration, hence the need for the new extended definition. Figure S14. Sliding energy curves of the antiparallel stacked homogeneous MoS2 bilayer calculated for vertically flexible shifts along the armchair direction using density functional theory (full curves) and the approach (dashed curves) using the previous definition. 23 The origin of the DFT energy scale is set to the total energy of the ′ ( ′ ) stacking mode and the curve is normalized by the energy of the 2 (Δ = 2 − ′ ) stacking mode.